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HIGHLIGHTS 


►  Effective  thermal  conductivity  increased  with  increasing  compression  and  GDL  thickness. 

►  Core  thermal  conductivities  are  higher  than  bulk  values  by  a  factor  of  2.87—7.79. 

►  Compared  to  PTFE  addition,  high  porosity  surfaces  dominate  thermal  conductivity. 
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In  this  work,  the  effective  through-plane  thermal  conductivities  for  compressed,  anisotropic  and 
heterogeneous  polymer  electrolyte  membrane  (PEM)  fuel  cell  gas  diffusion  layers  (GDLs)  were  deter¬ 
mined  analytically  from  representative  physical  GDL  models,  which  were  informed  by  microscale 
computed  tomography  imaging  of  four  commercially  available  GDL  materials.  The  number  of  fibre-to- 
fibre  contact  points  and  corresponding  contact  areas  were  extracted  from  these  physical  models  as 
inputs  to  a  thermal  resistance  model.  It  was  found  that  the  effective  thermal  conductivity  increased  with 
increasing  GDL  thickness  (with  bulk  porosity  remaining  almost  constant).  The  analytical  model  was 
employed  to  determine  the  bulk  thermal  conductivity  as  well  as  the  thermal  conductivity  of  the  core 
region  of  the  material.  By  isolating  the  core  region  from  the  bulk,  a  better  understanding  of  the  effect  of 
the  heterogeneous  porosity  profiles  on  the  through-plane  thermal  conductivity  was  determined  and 
discussed.  Unlike  the  bulk  thermal  conductivity,  the  thermal  conductivity  of  the  core  region  was  not 
dependent  upon  the  material  thickness.  It  was  also  found  that  the  surface  transition  regions  of  the 
porosity  distributions  have  a  dominating  effect  over  the  addition  of  PTFE  in  impacting  the  overall  thermal 
conductivity. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Achieving  proper  heat  management  within  a  polymer  electro¬ 
lyte  membrane  (PEM)  fuel  cell  is  critical  for  improving  its  perfor¬ 
mance  and  lifetime.  Heat  produced  from  the  electrochemical 
reaction  and  water  phase  change  results  in  temperature  gradients 
across  a  single  cell  and  fuel  cell  stack  [1-5].  Careful  control  of  the 
temperature  throughout  the  cell  is  important  for  avoiding 
membrane  dehydration  and  degradation  [4].  The  temperature 
within  the  fuel  cell  affects  the  relative  humidity,  membrane  water 
content,  saturation  pressure,  and  reaction  kinetics,  which  are 
coupled  together  in  their  impact  on  the  overall  performance  of  the 
fuel  cell.  Due  to  the  coupled  nature  of  these  parameters,  numerical 
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modelling  is  a  powerful  means  of  designing  the  PEMFC,  and  in 
particular,  designing  the  microscale  features  of  the  individual 
components.  However,  accurately  modelling  the  heat  transfer  rates 
and  temperature  distributions  within  a  fuel  cell  requires  knowl¬ 
edge  of  the  thermal  transport  properties  of  individual  materials,  in 
particular,  the  effective  thermal  conductivity. 

The  main  path  for  heat  removal  from  the  PEM  fuel  cell 
membrane  to  the  current  collectors  is  through  the  gas  diffusion 
layer  (GDL),  otherwise  known  as  the  diffusion  medium  (DM)  or 
porous  transport  layer  (PTL),  and  the  rate  of  heat  removal  is 
therefore  largely  dependent  upon  the  thermal  transport  properties 
of  the  GDL.  Although  a  number  of  analytical  correlations  exist  for 
the  effective  thermal  conductivity  of  a  composite  material  with 
varying  geometries  [6],  the  GDL  thermal  conductivity  can  be  esti¬ 
mated  but  not  accurately  represented  by  a  single  analytical  corre¬ 
lation  due  to  the  anisotropic  [7]  and  heterogeneous  nature  [8]  of 
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the  fibres  within  the  GDL.  Through  a  compact  analytical  model  of 
the  GDL,  Sadeghi  et  al.  [9]  calculated  the  through-plane  effective 
thermal  conductivity  of  the  GDL  using  an  idealized  repeating  unit 
cell  to  represent  the  GDL  structure.  The  analytical  model  accounted 
for:  the  conduction  through  the  solid  fibrous  material  and  gas 
phase,  the  thermal  constriction  resistance  at  the  regions  of  fibre-to- 
fibre  contact,  and  gas  rarefaction  in  the  microgaps  to  determine  the 
effective  thermal  conductivity.  The  results  from  their  model  [9]  and 
parametric  study  showed  that  the  thermal  constriction  resistance 
dominates  the  total  thermal  resistance  in  the  determination  of  the 
effective  thermal  conductivity. 

The  work  of  Sadeghi  et  al.  [9]  was  extended  by  the  same  authors 
to  include  an  experimental  and  analytical  study  on  the  effective 
thermal  conductivity  and  thermal  contact  resistance  of  the  GDL 
under  a  compressive  load  [10].  Increases  in  the  bipolar  plate  pres¬ 
sure  were  found  to  increase  the  thermal  conductivity.  Further  to 
this  approach,  Sadeghi  et  al.  [11]  recently  reported  the  first 
combined  experimental  and  analytical  study  on  the  in-plane 
effective  thermal  conductivity  and  thermal  contact  resistance  of 
the  GDL.  Their  analytical  model  accounted  for  heat  conduction 
through  randomly  oriented  fibres,  the  contact  area  between  fibres, 
and  polytetrafluoroethylene  (PTFE)  covered  regions.  The  results 
from  their  analytical  model  and  experimental  work  indicated  that 
PTFE  has  little  impact  on  the  in-plane  thermal  conductivity.  The 
results  obtained  for  the  in-plane  thermal  conductivity  were  12 
times  higher  than  the  through-plane  thermal  conductivity  from 
their  previous  work. 

The  effective  thermal  conductivity  of  the  GDL  can  be  deter¬ 
mined  in-situ  from  the  temperature  of  an  operating  fuel  cell.  Ex- 
situ  experiments  for  the  GDL  are  more  common  than  in-situ 
experiments  due  to  the  challenges  associated  with  the 
complexity  of  coupled  processes  within  an  operating  fuel  cell  [12]. 
The  first  experimentally  determined  GDL  thermal  conductivity  was 
reported  by  Vie  and  Kjelstrup  in  2003  [13].  Khandelwal  and  Mench 
[14]  measured  the  through-plane  thermal  conductivity  of  a  dry 
Nation  membrane  and  Sigracet  and  Toray  carbon  paper  GDLs  in 
a  steady-state  ex-situ  experiment.  The  results  presented  by  the 
authors  revealed  that  the  thermal  conductivity  is  a  function  of  the 
GDL  PTFE  content.  The  thermal  conductivity  of  an  untreated 
Sigracet  GDL  was  two  times  larger  than  the  thermal  conductivity  of 
a  Sigracet  GDL  treated  with  20%  PTFE.  Burheim  et  al.  [12]  reported 
experimental  measurements  for  the  anisotropic  through-plane 
effective  thermal  conductivity  of  dry  and  wet  Nation  membranes 
and  compressed  GDLs  from  ex-situ  experiments.  Compression  was 
shown  to  cause  a  nearly  linear  increase  on  the  effective  thermal 
conductivity  with  applied  compression  load.  Burheim  et  al.  [15] 
recently  reported  an  extensive  study  on  the  effective  through- 
plane  thermal  conductivity  of  a  number  of  GDL  materials.  The 
authors  [15]  studied  the  effect  of  residual  water,  PTFE  content,  and 
compression  on  the  through-plane  thermal  conductivity  and  the 
thermal  contact  resistance  between  the  GDL  and  aluminium  bi¬ 
polar  plate.  The  through-plane  thermal  conductivity  was  found  to 
decrease  with  increasing  PTFE  content. 

Sadeghi  et  al.  [9]  provided  important  insight  into  the  through- 
plane  thermal  conductivity  of  the  GDL  with  an  analytical  model 
considering  the  case  of  an  untreated  and  periodically  ordered  GDL 
structure.  In  practice,  the  GDL  is  anisotropic  and  carbon  paper  GDL 
materials  are  composed  of  randomly  oriented  fibres.  Recent  work 
by  Fishman  et  al.  [8]  indicated  that  the  porosity  of  the  GDL  is 
heterogeneous.  Burheim  et  al.  [15]  attributed  experimentally  found 
trends  in  the  through-plane  thermal  conductivity  with  the  GDL 
thickness  of  carbon  paper  materials  to  this  heterogeneity.  A 
numerical  investigation  into  the  impact  of  GDL  porosity  heteroge¬ 
neity  on  effective  thermal  conductivity  would  be  highly  informa¬ 
tive  for  designing  GDL  materials  for  optimized  heat  transfer 


properties,  which  would  lead  to  overall  improved  fuel  cell 
performance. 

In  this  work,  an  analytical  model  is  used  to  determine  the 
effective  through-plane  thermal  conductivity  of  commercially 
available  GDL  materials  based  on  the  heterogeneous  porosity 
profiles  previously  published  in  [8,16].  The  effects  of  the  hetero¬ 
geneous  porosity  distribution,  GDL  compression,  and  PTFE  content 
on  the  effective  thermal  conductivity  are  investigated.  The  results 
of  the  model  are  compared  with  experimental  data  presented  by 
Burheim  et  al.  in  Ref.  [15]. 

2.  Model  development 

Our  approach  to  determining  the  effective  through-plane 
thermal  conductivity  of  the  GDL  consists  of  two  steps:  the  repre¬ 
sentative  physical  GDL  model  and  the  thermal  resistance  model. 
Based  on  the  findings  of  Sadeghi  et  al.  [9,10]  where  heat  transfer 
was  found  to  be  dominated  by  fibre-to-fibre  contact,  our  analytical 
model  will  only  consider  heat  transfer  through  fibre-to-fibre 
contact  points.  Using  inputs  of  porosity,  GDL  compression  pres¬ 
sure,  and  in  the  case  of  treated  materials,  PTFE  distributions,  the 
representative  physical  GDL  model  is  used  to  determine  the 
number  of  fibre-to-fibre  contact  points  throughout  the  domain  and 
the  contact  area  dimensions  of  the  contact  points.  These  dimen¬ 
sions  are  then  used  as  inputs  to  the  thermal  resistance  model  for 
determining  the  effective  thermal  conductivity  of  the  GDL. 


2.1.  Fibre  representation 


To  model  the  effective  through-plane  thermal  conductivity  of 
the  GDL,  a  representative  modelling  domain,  shown  in  Fig.  1,  with 
width,  Wt  and  length,  L,  was  constructed  for  a  number  of  randomly 
oriented  fibres,  Nt,  in  the  x—y  plane.  The  fibres  are  immersed  in 
quiesent  air,  each  with  a  length  and  diameter  of  l  and  d,  respec¬ 
tively.  Layers  of  the  fibres  are  stacked  vertically  through  the  entire 
thickness,  H,  of  the  GDL  modelling  domain.  The  geometric 
parameters  of  the  modelling  domain  are  shown  in  Table  1. 

The  through-plane  porosity  distributions  reported  in  Ref.  [8]  of 
four  GDL  materials  are  employed  in  this  work  to  establish  the 
number  of  fibres,  in  each  layer,  t,  of  the  modelling  domain.  The 
porosity  of  each  layer  can  be  found  from: 


i  _  ^solid 
Vtotal 


(1) 


where  Vsoiid  is  the  total  volume  of  the  solid  carbon  fibres,  and  Vtotai 
is  the  total  volume  of  modelling  domain.  Based  on  the  volume  of 
the  modelling  domain  used  and  the  prescribed  porosity  of  each 
layer,  the  number  of  fibres  required  in  each  layer,  Nt  is  determined 
from: 


Nt  =  (1 


,LWt 

iS 

4 


(2) 


where  d  and  l  are  the  diameter  and  length  of  the  carbon  fibre, 
respectively. 

In  our  model,  the  number  of  contact  points,  Cf,  for  a  single  fibre 
in  layer,  t,  and  the  adjacent  layer,  t+  1,  is  a  function  of  the  porosity 
of  both  layers  (i.e.  number  of  fibres  in  each  layer,  Nt )  and  the 
orientation  angle  of  the  fibre,  6.  The  fibres  are  oriented  in  the 
modelling  domain  with  a  uniform  distribution  of  randomly  chosen 
angles  between  0  and  90°  (in  the  X—Y  plane).  The  maximum 
number  of  contact  points  for  a  given  fibre  in  each  layer  occurs  when 
the  fibre  orientation  angle  is  90°,  as  shown  in  the  periodically 
ordered  fibre  arrangement  (Fig.  2).  The  minimum  number  of 
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Fig.l.  Schematic  illustrating  a  section  of  the  representative  physical  model,  which  consists  of  solid,  continuous  fibres  oriented  in  a  planar  direction  a)  X—Y  plane  orientation,  and,  b) 
isometric  view. 


contact  points  is  seen  when  6  is  0°  between  the  two  layers.  A  non¬ 
linear  relationship  was  determined  for  the  number  of  contact 
points  for  a  single  fibre  in  a  layer  with  porosity,  et,  based  on  its 
orientation  angle  (Fig.  3a).  The  average  number  of  contact  points  for 
a  single  fibre,  Cf,  in  a  layer  with  porosity,  was  determined  from 
this  relationship  and  is  used  as  an  input  into  the  model.  The  average 
number  of  contact  points,  Cfiavg,  was  found  to  vary  linearly  with 
porosity  and  is  interpolated  for  each  porosity  value,  et  (Fig.  3b), 
when  generating  the  thermal  resistance  model. 

2.2.  Compression 


smooth  cylinders  is  used  to  define  the  elliptical  contact  region, 
described  by  the  major  and  minor  semi  axis,  a  and  b.  The  detailed 
formulation  for  this  theory  can  be  found  in  Sadeghi  et  al.  [9]  and 
Johnson  [17]. 

It  is  assumed  that  fibres  will  not  bend  and  new  contact  points 
will  not  be  created  under  compression.  Flowever,  under  compres¬ 
sion,  it  is  assumed  that  the  existing  contact  area  between  two  fibres 
will  increase.  The  contact  load,  F,  for  each  contact  point  can  be 
expressed  in  terms  of  the  pressure  within  the  GDL,  Pgdl.  the  cross- 
sectional  area  of  the  unit  cell  defined  by  L  and  W,  and  the  number  of 
contact  points  for  a  given  layer,  Ct  [9]: 


In  practice  the  PEM  fuel  cell  is  compressed  during  assembly.  It  is 
important  to  consider  the  effect  of  this  compression  on  both  the 
individual  carbon  fibres  and  on  the  overall  GDL.  Flere,  it  assumed 
that  the  fibres  do  not  bend  under  compression,  and  the  only 
deformation  of  the  fibres  is  at  the  contact  area  between  two  fibres. 
The  application  of  a  load  will  create  a  contact  area  between  two 
fibres  that  is  finite  but  small  when  compared  with  the  dimensions 
of  the  fibres  [17].  The  overall  thickness  of  the  GDL  will  change  due 
to  compression.  These  considerations  are  described  in  detail  in  the 
following  subsections. 


Fmax,t 


PgdlLW 

Ct 


(3) 


2.2.2.  GDL  compression 

In  this  work,  we  employed  the  porosity  profiles  and  material 
thicknesses  for  the  uncompressed  materials  measured  in  our 
group’s  previous  work  [8].  To  account  for  the  compression  of  the 
GDL,  we  employed  Hooke’s  law  to  determine  the  compressed 
thickness  of  the  material.  Assuming  the  deformation  of  the  GDL  is 
elastic,  ?,  the  strain  from  the  GDL  compression  is  given  by: 


2.2.1.  Fibre  compression 

Based  on  the  assumption  that  both  fibres  are  smooth,  a  smooth, 
non-conforming  contact  area  between  the  two  fibres  is  formed, 
which  is  a  function  of  the  force  applied  on  the  fibres,  F,  the  Young’s 
modulus  of  the  fibre  (210  GPa),  Ff,  and  the  angle  between  the  two 
fibres,  6.  The  material  properties  of  the  fibres  are  shown  in  Table  2. 

When  cylindrical  fibres  contact  each  other  eccentrically,  the 
contact  region  is  considered  elliptical  [17].  The  Hertzian  theory  of 
contact  is  used  to  predict  the  shape  of  contact  between  solids  and 
how  it  grows  under  an  increasing  load  [17].  Following  the  analytical 
modelling  approach  for  the  GDL  presented  by  Sadeghi  et  al.  [9],  the 
application  of  the  Hertzian  contact  theory  for  non-conforming 


Table  1 

Unit  cell  geometry  properties  employed  to  generate  the  representative  physical 
model. 


Fibre  diameter,  d 

Fibre  length,  l 

Unit  cell  width,  W 

Unit  cell  length,  L 

7.32  pm 

325  pm 

1500  pm 

1625  pm 

?  =  ^gdl/^gdl  (4) 

where  Fgdl  is  the  Young’s  modulus  of  the  GDL  material.  The 
compressed  thickness,  Hc,  is  determined  from: 


Fig.  2.  Schematic  illustrating  a  periodically  ordered  fibre  arrangement  with  6  =  90°. 
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Fig.  3.  (a)  Contact  points  per  fibre  in  layer,  t,  within  the  defined  unit  cell  with  a  porosity  of  80%  as  a  function  of  the  orientation  angle,  and,  (b)  average  number  of  contact  points  as 
a  function  of  porosity. 


where  H0  is  the  uncompressed  thickness  of  the  GDL.  A  value  of 
17.9  MPa  for  the  Young’s  modulus  of  Toray  carbon  paper  was  used 
[18]. 

The  analytical  model  assumes  that  compression  is  applied 
uniformly  within  the  GDL.  It  is  important  to  note  that,  this  is  not  the 
case  during  fuel  cell  operation  as  there  will  be  higher  compression 
in  areas  under  the  lands  of  the  bipolar  plate  than  in  areas  under  the 
channels  [12]. 

2.3.  PTFE  treatment 


through-plane  direction  of  the  domain  is  determined  for  this 
investigation  from  supplier  specifications  combined  with  the 
measurements  reported  in  Ref.  [20]. 

To  fully  define  the  modelling  domain,  the  number  of  fibres,  Nt ,  as 
well  as  the  weight  fraction  of  PTFE,  At,  must  be  prescribed  at  each 
through-plane  position  (layer)  of  the  GDL.  The  through-plane 
porosity  distributions  reported  in  Refs.  [8,16]  were  employed; 
however,  Fishman  et  al.  reported  porosity  distributions  that  accoun¬ 
ted  for  both  PTFE  and  carbon  fibres.  Therefore,  in  order  to  extract  Nt,  it 
is  first  necessary  to  determine  the  porosity  at  each  through-plane 
position  of  the  material  in  the  absence  of  PTFE  addition  (fibre 
porosity),  ro ,t>  by  using  Eq.  (5).  Once  Nt  is  extracted,  the  same  proce¬ 
dure  for  Fibre  Representation  and  Fibre  Compression  outlined  in 
Sections  2.1  and  2.2  for  untreated  GDL  materials  can  be  followed. 


For  the  GDL  materials  considered  that  have  PTFE  treatment,  the 
method  of  incorporating  the  PTFE  in  our  analytical  model  is 
adopted  from  previous  modelling  of  the  effective  thermal 
conductivity  in  the  in-plane  direction  by  Sadeghi  et  al.  [11].  The 
bulk  porosity  with  the  addition  of  PTFE,  £,  is  calculated  with  the 
general  form  of  the  equation: 


£ 


£o  ~P 


A(1  -  to) 
1  -A 


(6) 


where  £o  is  the  bulk  porosity  of  the  material  before  PTFE  is  added,  A 
is  the  weight  fraction  of  PTFE,  and  p  is  the  density  ratio  between 
carbon  fibre  and  PTFE  [19].  A  value  of  0.9  is  used  for  p  [18]. 

In  recent  experimental  work  by  Rofaiel  et  al.  [20],  the  relative 
through-plane  distribution  of  PTFE  in  Toray-TGP-H-090  was 
quantified  using  scanning  electron  microscopy  (SEM)  energy 
dispersive  X-ray  spectrometry  (EDS)  imaging.  Results  from  this 
work  revealed  an  almost  symmetrical  presence  of  PTFE  in  the  paper 
GDL  material  with  a  lower  accumulation  in  the  centre  region  and 
peaks  of  accumulation  towards  the  surfaces.  The  experimentally 
measured  PTFE  distributions  from  the  previously  published  work  of 
Rofaiel  et  al.  [20]  are  employed  in  this  investigation  to  inform  the 
non-uniform  distribution  of  PTFE. 

The  total  weight  fraction  of  PTFE  in  a  GDL  sample,  A,  is  known 
a  priori  from  supplier  specifications,  and  when  combined  with  the 
relative  through-plane  distribution  of  PTFE  provided  by  Rofaiel 
et  al.  [20],  the  heterogeneous  distribution  of  PTFE  in  the  through- 
plane  direction  of  the  GDL  can  be  extracted.  In  other  words,  the 
non-uniform  weight  fraction  of  PTFE  at  each  layer,  At,  in  the 


2.4.  Thermal  model 

The  effective  thermal  conductivity  is  a  porous  media  property  that 
accounts  for  the  contributions  of  the  thermal  conductivity  of  each 
phase  present  [21].  For  the  determination  of  the  effective  thermal 
conductivity  of  the  GDL,  the  two  phases  that  are  generally  considered 
are  solid  carbon  fibres  and  gaseous  air  with  values  of  thermal 
conductivity  of  120  W/m  I<  and  0.03  W/m  K,  respectively  [22]. 

For  the  analytical  model  presented  here,  it  is  assumed  that  the 
only  heat  transfer  is  steady-state,  one-dimensional  (1-D)  conduc¬ 
tion  in  the  through-plane  direction  of  the  GDL.  With  the  calculation 
of  the  Grashof  and  Peclet  numbers,  Ramousse  [7]  showed  that 
natural  convection  and  convective  heat  transfer  are  negligible 
compared  with  conduction  in  the  GDL.  Radiative  heat  transfer  can 
be  neglected  for  temperatures  below  1000  K  [23],  which  is  well 
above  the  operating  range  of  a  PEM  fuel  cell.  Therefore,  heat  is 
transferred  through  the  modelling  domain  from  fibre-to-fibre 
contact  only.  The  thermal  resistance  in  the  conduction  along 
a  fibre  is  ignored,  as  it  is  assumed  to  be  negligible  compared  with 
the  thermal  constriction  and  spreading  resistances  [9]. 

A  thermal  resistance  network  or  circuit  can  be  constructed  for  1- 
D  heat  transfer  with  no  internal  energy  generation  and  with 
constant  properties  [24].  The  thermal  resistance,  Rt)C0 nd»  for 
conduction  in  a  plane  wall  is  defined  by: 


Kr.cond  -  q  ~  ]<A 


Table  2 

Carbon  fibre  properties  employed  for  representative  physical  model. 


Thermal  conductivity,  k 

Poisson’s  ratio,  v 

Young’s  modulus,  Ef 

120  W/m  I<  [7] 

0.3  [9] 

210  GPa  [9] 

where  AT  is  the  difference  in  temperature  across  the  wall,  q  is  the 
heat  flux,  L  is  the  length  of  the  plane  wall,  k  is  the  thermal  conduc¬ 
tivity  of  the  wall,  and  A  is  the  cross  sectional  area  of  the  wall.  An 
equivalent  thermal  circuit  with  thermal  resistances  in  parallel  and 
series  is  analogous  to  an  electrical  circuit  governed  by  Ohm’s  law. 
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Fig.  4.  (a)  Equivalent  thermal  resistance  network  of  a  single  layer  with  Ct  number  of 
contact  points,  and  (b)  the  equivalent  thermal  resistance  network  through  the  GDL 
thickness  with  A  number  of  layers. 


The  dominant  thermal  resistance  for  heat  conduction  through 
the  GDL  is  the  thermal  constriction  and  spreading  resistance  [10]. 
The  thermal  spreading  resistance,  RSPj  is  equivalent  to  the  thermal 
constriction  resistance,  Rco,  and  accounts  for  the  thermal  energy 
that  is  transferred  between  the  two  fibres  at  the  contact  interface. 
When  the  dimensions  of  the  contact  area  are  very  small  compared 
with  the  dimensions  of  the  contacting  bodies,  the  heat  transfer 
through  the  contact  area  is  constrained,  and  the  solution  can  be 
modelled  as  a  heat  source  on  a  half-space  [25].  Heat  entering  the 
half-space  is  constricted  to  flow  through  the  small  contact  area 
(thermal  constriction  resistance),  and  heat  leaving  the  half-space 
spreads  out  from  the  contact  area  (thermal  spreading  resistance 

[25] .  In  the  case  of  the  GDL,  the  contact  area  between  the  two  fibres 
is  considered  as  the  heat  source,  and  the  much  larger  carbon  fibres 
are  considered  as  the  half  space. 

The  thermal  constriction  resistance  is  then  a  function  of  the 
elliptical  contact  area  calculated  with  the  Hertizan  contact  theory 
and  described  by  the  major  and  minor  semi-axis  of  the  elliptical 
contact  area,  a  and  b ,  respectively.  Using  the  work  of  Yovanovich 

[26]  that  resulted  in  half-space  solutions  for  different  contact 


interfaces,  Sadeghi  et  al.  [9]  expressed  the  thermal  constriction 
resistance  analytically  to  be: 


^co  — 


27 rksa 


-m 


(8) 


where  K(rj)  is  the  complete  elliptical  integral  and  a  function  of  the 
contact  area  dimensions,  a  and  b.  The  total  thermal  resistance  at 
a  fibre  contact  point,  Rc p,  is  given  by: 


Rcp  =  [Rco  +  Rsp]-1  (9) 

The  total  thermal  resistance  for  each  layer,  Rt  (Fig.  4a),  within  the 
modelling  domain  can  be  found  from  the  summation  of  each 
individual  thermal  resistance  of  each  contact  point  in  parallel: 


Rt 


(10) 


The  total  resistance  in  the  modelling  domain,  Ktotai  (Fig.  4b),  can 
be  found  as  a  series  summation  of  the  thermal  resistances  of  each 
layer,  for  a  total  of  A  layers  in  the  domain: 

Ktotal  =  X>  HD 

0 


The  effective  thermal  conductivity  can  be  found  from  the  total 
thermal  resistance: 


/<eff 


^•total  _  H  t 
Rtota\A  ^totalWI- 


(12) 


When  considering  GDL  materials  treated  with  PTFE,  the 
approach  outlined  above  is  followed,  and  only  heat  transfer 
through  fibre-to-fibre  contact  is  considered.  The  PTFE  distribution 
within  the  GDL  can  vary  with  method  of  application,  but  for  the 
purpose  of  this  analytical  model  only  PTFE  that  has  accumulated 
at  fibre  contact  points  will  be  considered.  An  equivalent  thermal 
resistance  at  each  fibre-to-fibre  contact  point  is  constructed  to 
include  the  addition  of  PTFE.  The  location  of  PTFE  at  a  contact 
point  and  the  equivalent  thermal  resistance  diagram  is  given  in 
Fig.  5.  The  total  thermal  resistance  at  a  fibre  contact  point,  Rc p,  is 
given  by: 


Fig.  5.  (a)  Schematic  illustrating  a  single  fibre-to-fibre  contact  Point  for  treated  carbon  PaPer,  and  (b)  the  equivalent  thermal  resistance  network  for  (a). 
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Kcp  =  [(^PTFE)  1+(^co+^SP)  1+(^PTFE)  (13) 

The  thermal  resistance  of  the  PTFE  is  given  by: 

Rfi,e  = 

where  a  and  0  are  parameters  developed  by  Sadeghi  et  al.  to 
describe  PTFE  [11].  The  proposed  relations  from  Sadeghi  et  al.  [11] 
describe  the  geometry  of  the  portion  of  PTFE  that  has  accumulated 
at  a  contact  point  from  SEM  images  of  carbon  paper  GDL  with 
varying  PTFE  wt.  content.  The  thermal  conductivity  of  PTFE  is 
0.649  W/m  K  [11  ].  The  parameter,  a,  gives  the  radius  ratio  of  contact 
area  to  fibre.  A  value  of  0.1  was  estimated  from  SEM  images  by 
Sadeghi  et  al.  [11  ]  and  is  also  used  in  this  study.  The  parameter  ft  is 
the  radius  ratio  of  PTFE  to  carbon  fibre  and  varies  with  the  amount 
of  PTFE  with  the  following  relationship  in  the  through-plane 
direction  [11]: 

ft  =  0.25  +  3  (At  —  0.05)  (15) 

where  Xt  is  weight  fraction  of  PTFE  at  each  layer  in  the  through- 
plane  direction.  As  before,  the  total  thermal  resistance  for  each 
layer,  Rt ,  within  the  modelling  domain  can  be  found  from  the 
summation  of  each  individual  thermal  resistance  of  each  contact 
point  in  parallel.  The  total  resistance  in  the  modelling  domain,  ftotai, 
can  be  found  as  a  series  summation  of  the  thermal  resistances  of 
each  layer. 

3.  Results 

The  analytical  model  was  implemented  with  heterogeneous 
through-plane  porosity  distributions  for  Toray  carbon  paper  TGP- 
FI-030,  060,  090,  and  120  obtained  through  X-ray  microscale 
computed  tomography  (pCT)  experiments  conducted  in  Ref.  [8,16]. 
The  GDL  materials  considered  were  uncompressed  and  did  not 
have  a  microporous  layer  (MPL).  The  spatial  resolution  for  the  pCT 
data  gathered  is  2.44  pm  [8,16];  however,  for  the  purpose  of  this 
model,  measurements  for  porosity  are  interpolated  from  the 
experimentally  determined  data  set  at  every  7.32  pm  through  the 
thickness  of  the  GDL.  The  details  of  the  microscale  computed 
tomography  visualization  are  presented  in  Ref.  [8,16]. 

3.1.  Untreated  paper 

The  effective  through-plane  thermal  conductivity  obtained  from 
the  analytical  model  is  shown  in  Table  3  for  three  compression 
pressures,  Pbp.  The  effective  thermal  conductivity  results  for 
varying  compression  pressures  from  the  analytical  model  are  also 
shown  in  Fig.  6  and  compared  with  experimental  data  for  Toray 
carbon  paper  GDL  materials  presented  by  Burheim  et  al.  [15].  The 
results  from  the  analytical  model  are  in  agreement  within  11.6%, 
7.6%,  and  4.0%  of  the  experimentally  obtained  results  from  Ref.  [15] 
for  Toray  TGP-H-060, 090, 120,  respectively  averaged  over  the  range 
of  compression  pressures.  Similar  trends  can  be  seen  between  both 
sets  of  data,  as  described  below. 

The  effective  thermal  conductivity  is  observed  to  increase  with 
increasing  compression  pressure.  As  noted  by  Burheim  et  al.  [15], 
and  shown  in  Fig.  6,  the  effective  thermal  conductivity  increases 
almost  linearly  with  increasing  compression  pressure.  The  results 
from  the  analytical  model  also  display  another  trend  noted  in  the 
experimental  work  by  Burheim  et  al.  [15]  for  Toray  carbon  GDL 
material  with  varying  thickness.  As  shown  in  Fig.  6,  the  effective 
thermal  conductivity  increases  with  increasing  GDL  thickness 


Table  3 

Effective  thermal  conductivity  of  untreated  Toray  paper  GDLs  for  different 
compression  pressures. 


PcDL(kPa) 

H([i  m) 

£bulk,avg  (%) 

keif  (W/m  K) 

Toray  TGP-H-030 

460 

114.11 

82.50 

0.2165 

930 

111.03 

82.01 

0.3233 

1390 

108.02 

81.51 

0.4099 

Toray  TGP-H-060 

460 

213.96 

81.62 

0.4505 

930 

208.18 

81.11 

0.6039 

1390 

202.54 

80.58 

0.7327 

Toray  TGP-H-090 

460 

290.03 

82.00 

0.4268 

930 

282.2 

81.50 

0.5945 

1390 

274.55 

80.98 

0.7317 

Toray  TGP-H-120 

460 

349.46 

78.13 

0.6111 

930 

340.03 

77.52 

0.8329 

1390 

330.81 

76.89 

0.9851 

(even  as  the  average  bulk  porosity  remains  approximately 
constant)  for  all  three  GDL  materials  presented  by  Burheim  et  al. 
[15].  The  results  presented  in  Fig.  6  for  the  analytical  model  follow 
this  trend  except  for  Toray  TGP-FI-060  and  090.  The  effective 
thermal  conductivity  of  Toray  TGP-H-060  is  an  average  of  3%  higher 
than  Toray  TGP-H-090.  To  further  explore  the  reason  for  this  trend, 
the  effective  thermal  conductivity  for  each  GDL  material  is  pre¬ 
sented  as  a  function  of  through-plane  position  in  Fig.  7. 

As  shown  in  Fig.  7,  the  thermal  conductivity  through  the 
thickness  of  the  GDL  is  strongly  dependent  upon  the  local  value  of 
the  porosity,  et.  For  all  four  GDL  materials  presented,  the  Toray 
carbon  paper  materials  display  local  porosity  minima  and  maxima 
throughout  the  thickness  [8].  At  the  location  of  a  porosity 
minimum,  a  maximum  local  value  of  thermal  conductivity  is 
observed.  The  four  materials  investigated  have  a  local  maximum 
value  of  thermal  conductivity  between  2.5  and  3.0  W/m  K. 

Fishman  et  al.  [8]  noted  that  the  heterogeneous  porosity 
distributions  for  all  four  GDL  materials  are  distinct  and  consist  of 
three  segments:  two  transitional  surface  regions  and  a  core  region. 
The  transitional  surface  region  extends  linearly  between  the  outer 


Fig.  6.  Comparison  of  measured  average  effective  thermal  conductivity  with 
compression  pressure  and  with  values  from  Burheim  et  al.  [15]. 
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Fig.  7.  Through-plane  effective  thermal  conductivity  and  porosity  distributions  for  untreated  (a)  Toray  TGP-H-030,  (b)  Toray  TGP-H-060,  (c)  Toray  TGP-H-090,  and  (d)  Toray  TGP-H- 
120. 


surfaces  and  the  local  porosity  minima  and  the  core  region  is 
between  the  two  transitional  surface  regions  (see  Fig.  8)  [8].  After 
isolating  the  core  region  in  the  analytical  model,  the  effective 
thermal  conductivity  of  this  region  was  determined.  Core  region 
thermal  conductivity  results  are  presented  in  Table  4  for  a  single 
compression  pressure,  Pbp,  of  460  kPa  and  compared  with  the 
overall  bulk  effective  thermal  conductivities.  The  core  thermal 
conductivity  is  7.79,  2.87,  3.74,  and  3.31  times  higher  than  the 
overall  bulk  thermal  conductivity  for  Toray  TGP-H-030,  060,  090, 
and  120,  respectively.  Unlike  the  overall  bulk  thermal  conductivity, 
the  thermal  conductivity  of  the  core  region  is  not  dependent  upon 
the  material  thickness. 

3.2.  PTFE  treated  paper 

The  analytical  model  was  implemented  with  heterogeneous 
through-plane  porosity  distributions  for  Toray  carbon  paper 


Fig.  8.  Through-plane  effective  thermal  conductivity  and  porosity  distributions  for 
Toray  TGP-H-090  depicting  core  region  defined  by  Fishman  et  al.  [8]. 


TGP-H-060  obtained  Ref.  in  [8,16]  for  four  PTFE  weight  contents;  0, 
5,  10,  and  20  wt.%.  The  effective  through-plane  thermal  conduc¬ 
tivity  obtained  from  the  analytical  model  is  shown  in  Fig.  9a  at 
a  single  compression  pressure,  Pbp,  of  460  kPa.  The  results  in  Fig.  9a 
show  an  increase  in  the  effective  thermal  conductivity  initially  as 
PTFE  is  added  to  the  material  (from  0  to  5  wt.%)  and  then  decrease 
as  the  amount  of  PTFE  increases  above  5  wt.%.  The  results  displayed 
in  Fig.  9a  are  counter-intuitive,  as  it  was  expected  that  the  thermal 
conductivity  would  increase  with  increasing  PTFE.  As  PTFE  is 
added,  additional  pathways  in  the  through-plane  direction  are 
created  for  heat  transfer.  To  better  understand  the  decrease  in  the 
thermal  conductivity  with  increasing  PTFE,  an  alternate  method  of 
applying  PTFE  in  the  analytical  model  was  investigated. 

In  this  alternate  method,  a  single  through-plane  porosity 
distribution  of  Toray  TGP-060  with  0  wt.%  PTFE  was  employed.  The 
desired  amount  of  PTFE  was  subsequently  added  to  the  material 
following  the  approach  outlined  in  Section  2.4  at  each  layer  in  the 
modelling  domain.  The  effective  through-plane  thermal  conduc¬ 
tivity  obtained  from  the  analytical  model  is  shown  in  Fig.  9b  for  the 
alternate  method  at  a  single  compression  pressure,  Pbp,  of  460  kPa. 
The  results  in  Fig.  9b  show  that  the  effective  thermal  conductivity 


Table  4 


Effective  thermal  conductivity  of  untreated  Toray  paper  GDLs:  bulk  and  core  values. 


^bulk.avg  £core 

^eff,  bulk  (W/m  K) 

core  (W/m  K) 

Toray  TGP-H-030 

0.825  0.742 

0.2165 

1.686 

Toray  TGP-H-060 

0.816  0.794 

0.4505 

1.294 

Toray  TGP-H-090 

0.820  0.801 

0.4268 

1.596 

Toray  TGP-H-120 

0.781  0.760 

0.4519 

1.496 
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Fig.  9.  Through-plane  effective  thermal  conductivity  with  0,  5, 10,  and,  20  wt.%  PTFE  using  (a)  four  distinct  porosity  profiles  presented  in  Fishman  [16]  for  Toray  TGP-H-060  and  (b) 
a  single  porosity  profile  of  Toray  TGP-H-060  with  0  wt.%  PTFE  employed  to  generate  representative  physical  GDL  models  of  treated  GDL  materials. 


increases  with  increasing  PTFE  content.  An  increase  of  24.2%  with 
the  addition  of  20  wt.%  PTFE  is  observed.  These  results  will  be 
discussed  later  in  Section  4.2. 

4.  Discussion 

4.1.  Untreated  paper 

The  effective  thermal  conductivity  is  observed  to  increase  with 
increasing  compression  pressure  (Fig.  6).  This  result  has  been 
previously  shown  in  experimental  [10,12,15]  and  analytical  work 
[9,10]  and  can  be  explained  with  the  analytical  model  presented.  As 
the  compression  pressure  increases,  the  elliptical  contact  area 
between  two  fibres  increases.  Increases  in  the  contact  area  cause 
a  decrease  in  the  thermal  constriction  and  spreading  resistance,  and 
subsequently,  an  increase  in  the  overall  thermal  conductivity. 
Compression  will  also  cause  a  decrease  in  the  overall  thickness  of  the 
GDL  material  that  has  been  accounted  for  in  this  analytical  model.  It 
is  important  to  note  that  while  a  value  of  17.9  MPa  was  employed  in 
this  work  for  the  Young’s  modulus  of  the  GDL,  within  the  literature, 
there  is  a  large  variation  [18,27];  however,  the  model  was  not  found 
to  be  significantly  sensitive  to  the  compressed  thickness,  and  the 
trends  reported  in  this  work  were  not  affected  by  the  value  used. 

Fishman  et  al.  [8]  noted  that  the  heterogeneous  porosity 
distributions  for  all  four  GDL  materials  are  distinct  but  each  display 
three  distinct  segments:  two  transitional  surface  regions  and  a  core 
region.  For  the  thinnest  GDL  investigated,  Toray  TGP-H-030,  this 
transitional  surface  region  accounts  for  approximately  66%  of  its 
total  thickness  [8].  Fishman  et  al.  [8]  observed  that  the  transitional 
surface  region  accounts  for  45%,  33%,  and  28%  of  the  material 
thickness  for  Toray  TGP- FI-060,  090,  and  120,  respectively.  The 
variation  of  the  overall  effective  thermal  conductivity  with  GDL 
thickness  observed  in  Table  3  can  be  attributed  to  the  relative 
amount  of  transitional  surface  region.  The  thermal  conductivity 
appears  to  be  strongly  affected  by  the  higher  porosity  values  in  the 
transitional  surface  regions,  regardless  of  the  overall  bulk  porosity 
value  or  the  local  maximum  thermal  conductivity  value.  This  can  be 
further  explained  by  the  results  in  Table  4  that  compare  the  overall 
bulk  thermal  conductivity  (/<eff,buik)  with  the  thermal  conductivity 
of  the  core  region  (/<eff,core)-  The  thermal  conductivity  of  the  core 
region  is  not  dependent  upon  the  material  thickness  but  appears  to 
depend  on  the  porosity  of  the  material. 

The  discrepancy  between  our  analytical  results  and  the  exper¬ 
imental  results  from  literature  in  Fig.  6  stems  from  the  small 
variability  in  our  experimentally  obtained  porosity  profiles  of 
commercially  available  materials.  Batch-specific  variations 


associated  with  the  manufacturing  process  have  been  observed  in 
these  commercially  available  materials  [28].  In  this  work,  one 
sample  of  each  GDL  was  employed  as  an  input  in  our  analytical 
model.  Flere,  we  see  that  the  results  in  Fig.  7  from  the  analytical 
model  are  quite  sensitive  to  the  porosity  values  in  the  transitional 
regions,  or  locations  of  porosity  maximums. 

4.2.  PTFE  treated  paper 

The  addition  of  PTFE  in  the  GDL  was  expected  to  increase  the 
effective  through-plane  thermal  conductivity  for  a  single  compres¬ 
sion  value.  Even  though  the  thermal  conductivity  of  PTFE  is  signifi¬ 
cantly  lower  than  that  of  carbon  fibre,  the  PTFE  is  expected  to  provide 
an  additional  pathway  for  heat  transfer  at  a  fibre  contact  point  in  the 
analytical  model.  The  results  from  the  analytical  model  for  Toray  TGP- 
FI-060  (see  Fig.  9a)  show  an  initial  increase  in  the  thermal  conduc¬ 
tivity  with  the  addition  of  5  wt.%  PTFE,  followed  by  a  counter-intuitive 
decrease  with  increasing  PTFE  content  (from  5  to  20  wt.%). 

An  alternate  method  of  PTFE  application  in  the  analytical  model 
was  used  to  further  understand  this  trend.  This  method  uses 
a  single  heterogeneous  porosity  distribution  (0  wt.%  PTFE)  upon 
which  PTFE  is  added,  and  the  results  for  the  thermal  conductivity 
(Fig.  9b)  show  an  increase  in  the  thermal  conductivity  with 
increasing  PTFE  content  (0-20  wt.%).  It  was  observed  for  the 
untreated  paper  GDL  materials  (Fig.  6)  that  the  thermal  conduc¬ 
tivity  was  strongly  affected  by  the  heterogeneous  porosity  distri¬ 
butions  and  the  higher  porosity  values  in  the  transitional  surface 
regions,  regardless  of  the  overall  bulk  porosity  value  or  the  local 
maximum  thermal  conductivity  value.  The  porosity  profiles  of 
Toray  TGP-H-060  with  0,  5,  10,  and  20  wt.%  PTFE  presented  by 
Fishman  et  al.  [16]  are  shown  in  Fig.  10.  The  addition  of  PTFE  affects 
the  overall  shape  of  the  through-plane  porosity  distribution  of  the 
material,  the  slope  of  the  transitional  surface  regions,  and  the 
values  of  the  porosity  maxima  and  minima.  The  slopes  of  the 
transitional  surface  regions  increases  with  increasing  PTFE  wt. 
content  (Fig.  10).  It  is  expected  that  the  larger  gradients  in  the 
porosity  distribution  at  the  transitional  surface  regions,  with  their 
associated  high  porosities  and  high  thermal  resistances,  strongly 
impact  the  overall  thermal  conductivity  of  the  PTFE  treated  GDL. 
From  the  results  of  the  alternate  method  of  PTFE  application,  it  is 
observed  that  the  surface  transition  regions  of  the  porosity  distri¬ 
butions  dominate  over  the  addition  of  PTFE  in  their  impact  on  the 
overall  thermal  conductivity.  The  decrease  in  thermal  conductivity 
with  increasing  PTFE  content  (from  5  to  20  wt.%)  in  Fig.  9(a)  is 
therefore  attributed  to  the  dominating  thermal  resistances  of  the 
surface  transition  regions  of  the  GDL. 
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Fig.  10.  Through-plane  porosity  distributions  Toray  TGP-H-060  from  Fishman  et  al. 
[16]  employed  to  generate  representative  physical  GDL  models. 

5.  Conclusion 

An  analytical  model  to  determine  the  effective  thermal  conduc¬ 
tivity  of  the  GDL  is  presented.  Representative  physical  GDL  models 
informed  by  microscale  computed  tomography  imaging  of  four 
commercially  available  uncompressed  GDL  materials  [8,16]  were 
employed  to  define  the  thermal  model.  The  effect  of  the  heteroge¬ 
neous  porosity  distribution,  GDL  compression,  and  PTFE  content  on 
the  effective  thermal  conductivity  was  investigated.  The  model 
predictions  for  untreated  Toray  carbon  paper  GDL  materials  were 
compared  with  recent  experimental  work  by  Burheim  et  al.  [15] 
where  two  district  trends  were  noted  between  both  sets  of  data. 
The  first  trend  noted  was  an  almost  linear  increase  in  the  effective 
thermal  conductivity  with  increasing  bipolar  plate  compaction 
pressure.  The  effective  thermal  conductivity  was  also  seen  to 
increase  with  increasing  GDL  thickness  as  bulk  porosity  remained 
almost  constant.  This  trend  was  attributed  to  the  heterogeneous 
porosity  profiles  of  the  material  as  the  thermal  conductivity  appears 
to  be  strongly  affected  by  the  higher  porosity  values  in  the  transi¬ 
tional  surface  regions,  regardless  of  the  overall  bulk  porosity  value  or 
the  local  maximum  thermal  conductivity  value.  The  analytical 
modelling  approach  allows  this  trend  to  be  investigated  further  by 
isolating  the  core  region  and  studying  its  effective  thermal 
conductivity.  The  effective  thermal  conductivity  of  the  core  region 
was  found  to  be  independent  of  the  material  thickness. 

Two  methods  of  applying  PTFE  in  the  analytical  model  were 
investigated  to  isolate  the  impact  of  the  PTFE  addition  and  the 
through-plane  porosity  distribution  on  the  effective  thermal 
conductivity.  An  initial  increase  in  the  thermal  conductivity  with 
the  addition  of  5  wt.%  PTFE,  followed  by  a  decrease  with  increasing 
PTFE  content  (from  5  to  20  wt.%)  was  noted  when  four 


heterogeneous  porosity  distributions  were  employed  in  the 
analytical  model.  Flowever,  an  overall  increase  of  24.2%  in  the 
effective  thermal  conductivity  was  noted  with  the  addition  of 
20  wt.%  PTFE  when  alternate  method  of  PTFE  application  was 
employed,  thus  indicating  that  the  overall  thermal  conductivity  has 
a  strong  dependence  on  the  heterogeneous  porosity  distributions 
of  the  treated  GDL  materials  (in  particular,  the  surface  transition 
regions).  The  outcomes  of  this  work  provide  insight  into  domi¬ 
nating  effect  of  heterogeneity  and  anisotropy  of  the  GDL  on  the 
thermal  management  required  for  improved  PEMFC  performance. 
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